library(landscapemetrics)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.2.1     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.3
## ✓ tidyr   1.0.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.4.0
## ── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(ggrepel)
library(readxl)
library(stars)
## Loading required package: abind
## Loading required package: sf
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0

Восстанавливаем цветовую шкалу

typesdf = read_excel('data/types.xlsx')
types = typesdf$type

colormap = read_table2('data/colormap.clr', 
                       col_names = c('N', 'R', 'G', 'B', 'A', 'type')) %>% 
  filter(type %in% types)
## Parsed with column specification:
## cols(
##   N = col_double(),
##   R = col_double(),
##   G = col_double(),
##   B = col_double(),
##   A = col_double(),
##   type = col_double()
## )
pal = rgb(colormap$R/255, colormap$G/255, colormap$B/255, colormap$A/255)

Грузим данные

read_lc = function(file) {
  read_stars(file) %>% 
    select(type = 1) %>%
    mutate(typefac = factor(type, levels = types),
           urban = (type == 50),
           crops = (type == 40)) %>% 
    st_warp(crs = 32662)
}

salekhard = read_lc('data/SALEKHARD.tif')
moscow = read_lc('data/MOSCOW.tif')
rostov = read_lc('data/ROSTOV.tif')
ufa = read_lc('data/UFA.tif')
petro = read_lc('data/PETRO.tif')
surgut = read_lc('data/SURGUT.tif')
ustug = read_lc('data/USTUG.tif')
grozny = read_lc('data/GROZNY.tif')
voronezh = read_lc('data/VORONEZH.tif')

Визуализируем:

plot_lc = function(lc) {
  ggplot() +
    geom_stars(data = lc['typefac']) +
    scale_fill_manual(values = pal, breaks = types, labels = typesdf$class, drop = FALSE) +
    coord_sf(crs = st_crs(lc)) +
    xlab('LON') + ylab('LAT') +
    theme_bw() + theme(legend.position = 'bottom', legend.title=element_blank()) +
    guides(fill = guide_legend(ncol=4))
}

plot_lc(salekhard)
## Warning: Removed 4812 rows containing missing values (geom_raster).

plot_lc(moscow)
## Warning: Removed 4812 rows containing missing values (geom_raster).

plot_lc(rostov)
## Warning: Removed 4812 rows containing missing values (geom_raster).

plot_lc(ufa)
## Warning: Removed 4812 rows containing missing values (geom_raster).

plot_lc(petro)
## Warning: Removed 4812 rows containing missing values (geom_raster).

plot_lc(surgut)
## Warning: Removed 4812 rows containing missing values (geom_raster).

plot_lc(ustug)
## Warning: Removed 3208 rows containing missing values (geom_raster).

plot_lc(grozny)
## Warning: Removed 4812 rows containing missing values (geom_raster).

plot_lc(voronezh)
## Warning: Removed 4812 rows containing missing values (geom_raster).

Считаем метрики сложности:

lands = lst(САЛЕХАРД = salekhard, МОСКВА = moscow, `РОСТОВ-НА-ДОНУ` = rostov, УФА = ufa, ПЕТРОЗАВОДСК = petro,
            ШАРЬЯ = ustug, ГРОЗНЫЙ = grozny, ВОРОНЕЖ = voronezh, СУРГУТ = surgut)

(metrics = imap(lands, function(land, name) {
  tibble(name = name,
         ent = lsm_l_ent(land) %>% pull(value),
         condent = lsm_l_condent(land) %>% pull(value),
         joinent = lsm_l_joinent(land) %>% pull(value),
         mutinf = lsm_l_mutinf(land) %>% pull(value),
         relmutinf = mutinf / ent,
         urban = 100 * sum(land$urban, na.rm = TRUE) / length(land$urban),
         crops = 100 * sum(land$crops, na.rm = TRUE) / length(land$crops))
}) %>% bind_rows())
## # A tibble: 9 x 8
##   name             ent condent joinent mutinf relmutinf  urban   crops
##   <chr>          <dbl>   <dbl>   <dbl>  <dbl>     <dbl>  <dbl>   <dbl>
## 1 САЛЕХАРД        2.52   1.13     3.65  1.39      0.552  0.115  0     
## 2 МОСКВА          2.95   1.38     4.33  1.57      0.532 17.4   17.3   
## 3 РОСТОВ-НА-ДОНУ  1.49   0.561    2.05  0.927     0.623  3.20  72.4   
## 4 УФА             2.16   0.977    3.13  1.18      0.547  3.48  45.9   
## 5 ПЕТРОЗАВОДСК    2.35   0.945    3.30  1.41      0.599  0.196  0.985 
## 6 ШАРЬЯ           2.03   0.853    2.88  1.18      0.579  0.245  4.47  
## 7 ГРОЗНЫЙ         2.16   0.780    2.94  1.38      0.640  2.75  22.6   
## 8 ВОРОНЕЖ         1.86   0.759    2.62  1.10      0.592  3.09  65.5   
## 9 СУРГУТ          2.90   1.27     4.17  1.63      0.561  0.343  0.0543
ggplot(metrics, mapping = aes(ent, relmutinf, label = name)) +
  geom_point()

1 Пространственные данные

Читаем данные:

scales = c(200, 500, 1000)
regions = c('МОСКВА', 'ПЕТРОЗАВОДСК', 'САЛЕХАРД', 'УФА', 'РОСТОВ-ДОН')

lindata = lapply(scales, function(scale) { 
    read_delim(paste0('tables/resLines', scale, '.txt'), delim = ';', locale = locale(decimal_mark = ".")) %>% 
    select(-ncol(.)) %>% 
    mutate(Scale = scale,
           Region = stringr::str_sub(Region, 1, nchar(Region) - nchar(as.character(Scale))))
  }) %>% 
  bind_rows() %>% 
  filter(!stringr::str_detect(Name, 'clip')) %>% 
  mutate(Layer = stringr::str_sub(Name, 1, 3),
         Dim = stringr::str_sub(Name, 4, 6)) %>% 
  filter(!(Layer %in% c('rlh', 'rlf', 'veg'))) %>% 
  mutate_all(~replace(., is.nan(.), 0))
## Warning: Missing column names filled in: 'X17' [17]
## Parsed with column specification:
## cols(
##   Region = col_character(),
##   Name = col_character(),
##   PointsNumber = col_double(),
##   BendNumber = col_double(),
##   AveAngle = col_double(),
##   TotalArea = col_double(),
##   Length = col_double(),
##   AverageLength = col_double(),
##   AverageBaseLineLength = col_double(),
##   AverageHeight = col_double(),
##   AverageWidth = col_double(),
##   AverageArea = col_double(),
##   AtrbCount = col_double(),
##   UniqueValues = col_double(),
##   AveUniqueValues = col_double(),
##   NumberOfObj = col_double(),
##   X17 = col_logical()
## )
## Warning: Missing column names filled in: 'X17' [17]
## Parsed with column specification:
## cols(
##   Region = col_character(),
##   Name = col_character(),
##   PointsNumber = col_double(),
##   BendNumber = col_double(),
##   AveAngle = col_double(),
##   TotalArea = col_double(),
##   Length = col_double(),
##   AverageLength = col_double(),
##   AverageBaseLineLength = col_double(),
##   AverageHeight = col_double(),
##   AverageWidth = col_double(),
##   AverageArea = col_double(),
##   AtrbCount = col_double(),
##   UniqueValues = col_double(),
##   AveUniqueValues = col_double(),
##   NumberOfObj = col_double(),
##   X17 = col_logical()
## )
## Warning: Missing column names filled in: 'X17' [17]
## Parsed with column specification:
## cols(
##   Region = col_character(),
##   Name = col_character(),
##   PointsNumber = col_double(),
##   BendNumber = col_double(),
##   AveAngle = col_double(),
##   TotalArea = col_double(),
##   Length = col_double(),
##   AverageLength = col_double(),
##   AverageBaseLineLength = col_double(),
##   AverageHeight = col_double(),
##   AverageWidth = col_double(),
##   AverageArea = col_double(),
##   AtrbCount = col_double(),
##   UniqueValues = col_double(),
##   AveUniqueValues = col_double(),
##   NumberOfObj = col_double(),
##   X17 = col_logical()
## )
pntdata = lapply(scales, function(scale) { 
    read_delim(paste0('tables/resPoints', scale, '.txt'), delim = ';', locale = locale(decimal_mark = ".")) %>% 
    select(-ncol(.)) %>% 
    mutate(Scale = scale,
           Region = stringr::str_sub(Region, 1, nchar(Region) - nchar(as.character(Scale))))
  }) %>% 
  bind_rows() %>% 
  mutate(Layer = stringr::str_sub(Name, 1, 3),
         Dim = stringr::str_sub(Name, 4, 6)) %>% 
  filter(!(Layer %in% c('rlh', 'rlf', 'veg')))
## Warning: Missing column names filled in: 'X7' [7]
## Parsed with column specification:
## cols(
##   Region = col_character(),
##   Name = col_character(),
##   AtrbCount = col_double(),
##   UniqueValues = col_double(),
##   AveUniqueValues = col_double(),
##   NumberOfObj = col_double(),
##   X7 = col_logical()
## )
## Warning: Missing column names filled in: 'X7' [7]
## Parsed with column specification:
## cols(
##   Region = col_character(),
##   Name = col_character(),
##   AtrbCount = col_double(),
##   UniqueValues = col_double(),
##   AveUniqueValues = col_double(),
##   NumberOfObj = col_double(),
##   X7 = col_logical()
## )
## Warning: Missing column names filled in: 'X7' [7]
## Parsed with column specification:
## cols(
##   Region = col_character(),
##   Name = col_character(),
##   AtrbCount = col_double(),
##   UniqueValues = col_double(),
##   AveUniqueValues = col_double(),
##   NumberOfObj = col_double(),
##   X7 = col_logical()
## )
poldata = lapply(scales, function(scale) { 
    read_delim(paste0('tables/resPolygons', scale, '.txt'), delim = ';', locale = locale(decimal_mark = ".")) %>% 
    select(-ncol(.)) %>% 
    mutate(Scale = scale,
           Region = stringr::str_sub(Region, 1, nchar(Region) - nchar(as.character(Scale))))
  }) %>% 
  bind_rows() %>% 
  mutate(Layer = stringr::str_sub(Name, 1, 3),
         Dim = stringr::str_sub(Name, 4, 6)) %>% 
  filter(!(Layer %in% c('rlh', 'rlf', 'veg')))
## Warning: Missing column names filled in: 'X12' [12]
## Parsed with column specification:
## cols(
##   Region = col_character(),
##   Name = col_character(),
##   PointsNumber = col_double(),
##   TotalArea = col_double(),
##   TotalPerimetr = col_double(),
##   AveArea = col_double(),
##   AvePerimetr = col_double(),
##   AtrbCount = col_double(),
##   UniqueValues = col_double(),
##   AveUniqueValues = col_double(),
##   NumberOfObj = col_double(),
##   X12 = col_logical()
## )
## Warning: Missing column names filled in: 'X12' [12]
## Parsed with column specification:
## cols(
##   Region = col_character(),
##   Name = col_character(),
##   PointsNumber = col_double(),
##   TotalArea = col_double(),
##   TotalPerimetr = col_double(),
##   AveArea = col_double(),
##   AvePerimetr = col_double(),
##   AtrbCount = col_double(),
##   UniqueValues = col_double(),
##   AveUniqueValues = col_double(),
##   NumberOfObj = col_double(),
##   X12 = col_logical()
## )
## Warning: Missing column names filled in: 'X12' [12]
## Parsed with column specification:
## cols(
##   Region = col_character(),
##   Name = col_character(),
##   PointsNumber = col_double(),
##   TotalArea = col_double(),
##   TotalPerimetr = col_double(),
##   AveArea = col_double(),
##   AvePerimetr = col_double(),
##   AtrbCount = col_double(),
##   UniqueValues = col_double(),
##   AveUniqueValues = col_double(),
##   NumberOfObj = col_double(),
##   X12 = col_logical()
## )
nuniq = length(unique(poldata$Region))

inters = read_delim('tables/topLines.txt', delim = ';', locale = locale(decimal_mark = ".")) %>% 
  select(-ncol(.)) %>% 
  mutate(Scale = as.character(rep(scales, each = n()/3)),
         Region = stringr::str_sub(Region, 1, nchar(Region) - nchar(Scale)))
## Warning: Missing column names filled in: 'X3' [3]
## Parsed with column specification:
## cols(
##   Region = col_character(),
##   Intersections = col_double(),
##   X3 = col_logical()
## )

Подсчитаем суммарное количество социально-экономических точек:

linpts = lindata %>% 
  group_by(Region, Scale) %>% 
  summarise(linnpts = sum(PointsNumber))
polpts = poldata %>% 
  group_by(Region, Scale) %>% 
  summarise(polnpts = sum(PointsNumber),
            area = 1e-6 * sum(TotalArea * (Layer == 'adm')))
pntpts = pntdata %>% 
  group_by(Region, Scale) %>% 
  summarise(pntnpts = sum(NumberOfObj))


stats = bind_cols(linpts, polpts, pntpts) %>% 
  select(name = Region, Scale, pntnpts, linnpts, polnpts, area) %>% 
  mutate(total = pntnpts + linnpts + polnpts,
         density = total / area) %>% 
  ungroup() %>% 
  mutate(name = ifelse(name == 'РОСТОВ-ДОН', 'РОСТОВ-НА-ДОНУ',
                       ifelse(name == 'УСТЮГ', 'ШАРЬЯ', name)))

Соединяем

compare = left_join(stats, metrics)
## Joining, by = "name"

Корреляционный и регрессионный анализ

tab = compare %>% 
  filter(Scale == 1000) 

ggplot(tab, mapping = aes(urban, density, label = name)) +
  geom_point() +
  geom_smooth(method = 'lm') +
  geom_text_repel(
    segment.size  = 0.3,
    segment.color = "grey50",
    direction     = "y",
    hjust         = 0
  ) +
  scale_x_log10()

ggplot(tab, mapping = aes(crops, density, label = name)) +
  geom_point() +
  geom_smooth(method = 'lm') +
  geom_text_repel(
    segment.size  = 0.3,
    segment.color = "grey50",
    direction     = "y",
    hjust         = 0
  ) +
  scale_x_log10()
## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Removed 1 rows containing non-finite values (stat_smooth).

ggplot(tab, mapping = aes(crops + urban, density, label = name)) +
  geom_point() +
  geom_smooth(method = 'lm') +
  geom_text_repel(
    segment.size  = 0.3,
    segment.color = "grey50",
    direction     = "y",
    hjust         = 0
  ) +
  scale_x_log10()

cor.test(tab$density, log(tab$urban + 1))
## 
##  Pearson's product-moment correlation
## 
## data:  tab$density and log(tab$urban + 1)
## t = 4.9167, df = 7, p-value = 0.00172
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5214060 0.9746945
## sample estimates:
##       cor 
## 0.8805971
cor.test(tab$density, log(tab$crops + 1))
## 
##  Pearson's product-moment correlation
## 
## data:  tab$density and log(tab$crops + 1)
## t = 3.0497, df = 7, p-value = 0.01859
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1830774 0.9452803
## sample estimates:
##       cor 
## 0.7553583
cor.test(tab$density, log(tab$crops + 1) + log(tab$urban + 1))
## 
##  Pearson's product-moment correlation
## 
## data:  tab$density and log(tab$crops + 1) + log(tab$urban + 1)
## t = 4.3883, df = 7, p-value = 0.003202
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4458400 0.9692525
## sample estimates:
##       cor 
## 0.8563911
lm(density ~ log(urban + 1), tab)
## 
## Call:
## lm(formula = density ~ log(urban + 1), data = tab)
## 
## Coefficients:
##    (Intercept)  log(urban + 1)  
##         0.0646          0.1251
lm(density ~ log(crops + 1), tab)
## 
## Call:
## lm(formula = density ~ log(crops + 1), data = tab)
## 
## Coefficients:
##    (Intercept)  log(crops + 1)  
##        0.06255         0.05715
lm(density ~ log(crops + urban + 1), tab)
## 
## Call:
## lm(formula = density ~ log(crops + urban + 1), data = tab)
## 
## Coefficients:
##            (Intercept)  log(crops + urban + 1)  
##                0.04241                 0.06152

Зависимость плотности от сложности ландшафта:

ggplot(tab, mapping = aes(joinent, density, label = name)) +
  geom_point() +
  geom_text_repel(
    segment.size  = 0.3,
    segment.color = "grey50",
    direction     = "y",
    hjust         = 0
  )